variables

ts=(t0:stepsize:tend)'; %sampling times
N=length(ts); % Anzahl Schritte

testinp=@(t) sin(omega*t.*t); %input function
u=testinp(ts); %inputs at times ts

ddd=[1;2;3;4;3;2;7]*0.1;
%R=spdiags(repmat(ddd,ceil(N/length(ddd)),1),0,N,N); -> in variables
%ausgelagert
disturbance=R*random('norm',0,ones(size(u)));
Q=inv(R);

testt=tf([ucoeff;0].',[1;-ycoeff].',stepsize,'Variable','z');
% SYS = TF(NUM, DEN, TS, 'PropertyName1', 'PropertyValue1')
%  creates a discrete-time transfer function with
%  sample time TS (set TS=-1 if the sample time is undetermined), 
%  with numerator(s) NUM and denominator(s) DEN.  
%  The output SYS is a TF object.
%------------------
%testt: 
% 
%Transfer function:
%    0.73 z^3 + 3 z^2 + 7 z
%------------------------------
%z^3 - 0.15 z^2 - 0.27 z + 0.43
% 
%Sampling time: 0.01
%-----------------

y=dlsim(testt,u,disturbance,ts);
% dlsim(sys, u, v, t) returns the time response of the lti model 
% sys', where sys is of the form y(z) = B(z)/A(z) * u(z) und sys' 
% is y(z) = B(z)/A(z) + 1/A(z) * v(z), that is, a _d_isturbance v
% (hence dlsim) is added in every step. We need this specific 
% model structure because it is the model structure identified by
% the basic least squares model identification algorithm, that we
% have programmed and want to test.
% At the moment t is ignored.

phat = lsidentify(y,u,yl,ul,Q);

% Kontrolle: es kommt ca. coeff=[ycoeff,ucoeff] raus
phat
phat-coeff
norm(phat-coeff)

plot(ts,y,'m',ts,u,'g') %Input und Output des Systems, das wir schaetzen
